// ExpWaveProb.C

#include <iostream>
#include <fstream>

#include "Epetra_SerialComm.h"
#include "Epetra_Map.h"

#include "ExpWaveProb.H"

using namespace std;

ExpWaveProb::ExpWaveProb(int m, int n, double p0, double a, double b)
{
  M = m;
  N = n;
  P0 = p0;
  A = a;
  B = b;
  firstFill = true;

  Epetra_SerialComm comm;
  map = new Epetra_Map(M*N, 0, comm);
  initial = new Epetra_Vector(*map);
  J = new Epetra_CrsMatrix(Copy, *map, 9);

  double lambda = computeLambda();
  initQ(lambda);
  computeInitial(lambda);
}

ExpWaveProb::ExpWaveProb(char *solfile)
{
  firstFill = true;

  ifstream fin(solfile); 
  if (fin.fail()) { 
    cerr << "Failed to open initial solution file, exiting" << endl; 
    exit(1); 
  }
  char header[32];
  fin >> header; 
  if (strcmp(header, "exp")) { 
    cerr << "Error: initial solution file provided does not begin with \"exp\", exiting" << endl; 
    exit(1); 
  }
  fin >> M;
  fin >> N;
  fin >> P0;
  fin >> Q;
  fin >> A;
  fin >> B;

  Epetra_SerialComm comm;
  map = new Epetra_Map(M*N, 0, comm);
  initial = new Epetra_Vector(*map);
  J = new Epetra_CrsMatrix(Copy, *map, 9);

  double v;
  for (int i=0; fin >> v; i++)
    (*initial)[i] = v;
  fin.close();
}

ExpWaveProb::~ExpWaveProb()
{
  delete map;
  delete initial;
  delete J;
}

double ExpWaveProb::computeLambda()
{
  double hp = fabs(P0) / M;
  double y;
  double intGamma = 0.0;
  double gammaMin = 0.0;
  
  //  Special insert for gamma = 0 case. 
  if (B == 0) // gamma = const
  {
  	if (A == 0)
  	{
    	//  double WaveProb::newton(double guess)
    	y = 1;
    	double fy = y - g*tanh(fabs(P0)/sqrt(y));
    	if (fabs(fy) > 1e-25)
      	{
		double dfy = 1 + (g*fabs(P0)*0.5)*pow(y, -1.5)*pow(sech(fabs(P0)/sqrt(y)), 2);
		double d = -fy/dfy;
		while (fabs(d)/(fabs(y) + 1e-5) > 1e-8)
	  	{
	    	y += d;
	    	fy = y - g*tanh(fabs(P0)/sqrt(y));
	    	dfy = 1 + (g*fabs(P0)*0.5)*pow(y, -1.5)*pow(sech(fabs(P0)/sqrt(y)), 2);
	    	d = -fy/dfy;
	  	}
      }	
      cerr << endl << "lambda* = " << y << "..." << endl;
      return y;
  	}
   // Special insert for constant  gamma.  Can't seem to converge 
   // using the general gamma case for constant gamma, especially at extreme ends
   // The plots of f(y) tell the story;  we need to choose the initial guess 
   // accordingly.  Altered the dfy term since otherwise it seems hard to
   // converge onto a correct answer. Changed pow(9.8 + A*sqrt(y), 2) to
   // 9.8 +  A*y.  Comes VERY close and seems more robust.
  	else
    {
    	if (A < 0) y = 1e-2;
    	else y  = 100;
      	double fy = tanh((sqrt(y) - sqrt(y + 2*P0*A))/A) - y/(9.8+A*sqrt(y));
      	if (fabs(fy) > 1e-25)
		{
	  		double  dfy = pow(sech((sqrt(y)-sqrt(y+2*P0*A))/A), 2)*(1/sqrt(y) - 1/sqrt(y+2*P0*A))/(2*A) - (9.8 + 0.5*sqrt(y)*A)/(9.8 +  A*y);
	  		double d = -fy/dfy;
	  		while (fabs(d)/(fabs(y) + 1e-5) > 1e-8)
	    	{
	    		y += d;
	      		fy = tanh((sqrt(y) - sqrt(y + 2*P0*A))/A) - y/(9.8+A*sqrt(y));
	      		dfy = pow(sech((sqrt(y)-sqrt(y+2*P0*A))/A), 2)*(1/sqrt(y) - 1/sqrt(y+2*P0*A))/(2*A) - (9.8 + 0.5*sqrt(y)*A)/(9.8 + A*y);
	      		d = -fy/dfy;
	    	}
	 	}	
	    cerr << endl << "lambda* = " << y << "..." << endl;
      	return y;  	
    }
  }
  //General, non-constant exponential gamma.  
  //We know that lambda* \in [-2*Gamma_min, lambda0].  In most cases, a good
  //start is lambda0.  IN the case of the Tetsu vorticity, a better guess is
  //-2*Gamma_min = 0
  else
  {	      
      	for (int r=0; r<M; ++r)
	 	{
	  		intGamma -= (Gamma(r*hp)+Gamma((r+1)*hp))*hp/2;
	  		if (intGamma < gammaMin) { gammaMin = intGamma; }
    	}  
      	y =  newton(-2*gammaMin+1e-2);
      	cerr << endl << "lambda0 = " << y << endl;
      	return y/100;
  }
}

void ExpWaveProb::initQ(double lambda)
{
  double hp = fabs(P0) / M;
  double intGamma = 0.0;
  double next;
  Q = 0.0;
  if (B == 0)
  {
  	//if (A == 0) Q = lambda + 2*g*fabs(P0)/sqrt(lambda);
  	//else Q = lambda + (4*g*fabs(P0))/(sqrt(lambda) + sqrt(lambda + 2*P0*A));
  	Q = lambda + (4*g*fabs(P0))/(sqrt(lambda) + sqrt(lambda + 2*P0*A));
  }
  else
  {
    for (int r=0; r<M-1; ++r)
	 {
	  	intGamma -= (Gamma(r*hp)+Gamma((r+1)*hp))*hp*0.5;
	  	Q -= (hp/2)/sqrt(lambda + 2*intGamma);
	  	next = intGamma - (Gamma((r+1)*hp)+ Gamma((r+2)*hp))*hp*0.5;
	  	Q -= (hp/2)/sqrt(lambda + 2*next);
	  }
      int r = M-1;
      intGamma -= (Gamma(r*hp) + Gamma((r+1)*hp))*hp*0.5;
      Q -= hp/sqrt(lambda + 2*intGamma);
 
      Q = -2*g*Q + lambda;
   }
   //Test out the value of Q0 = Q(0)  
   // double Q0 = 0.0;
   //	 intGamma = 0.0;
   //  for (int r=0; r<M-1; ++r)
   // {
   //	  	intGamma -= (Gamma(r*hp)+Gamma((r+1)*hp))*hp*0.5;
   //	  	Q0 = Q0- (hp/2)/sqrt(2*intGamma);
   //	  	next = intGamma - (Gamma((r+1)*hp)+Gamma((r+2)*hp))*hp*0.5;
   //	  	Q0 = Q0- (hp/2)/sqrt(2*next);
   //	  }
   //   int r = M-1;
   //   intGamma -= (Gamma(r*hp) + Gamma((r+1)*hp))*hp*0.5;
   //   Q0 = Q0- hp/sqrt(2*intGamma);
   //	  Q0 = -2*g*Q0;
   //  cerr << endl << "Q(0) = " << Q0 << "..." << endl;
	  
 
  cerr << endl << "Starting search for Q* at Q = " << Q << "..." << endl;
}


void ExpWaveProb::computeInitial(double lambda)
{
  double hp = fabs(P0)/M;

  if (B == 0) // gamma = 0
   {
   	  if (A == 0)
   		{
        	for (int r=0; r<M; ++r)
				for (int c=0; c<N; ++c)
	  				(*initial)[r*N + c] = (-hp*r - P0) / sqrt(lambda);
    	}
       else
    	{
      		for (int r=0; r<M; ++r)
				for (int c=0; c<N; ++c)
	  				(*initial)[r*N + c] = ( sqrt(lambda - 2*hp*r*A) - sqrt(lambda +  2*P0*A ))/A;
    	}
    }
    else // general gamma
  	{
    	double t2 = (Q - lambda)/(2*g);
    	double integral = 0;
    	double intGamma = 0;
    	double next;
    	for (int r=0; r<M-1; ++r)
    	{
      		intGamma -= (Gamma(r*hp)+ Gamma((r+1)*hp))*hp/2;
      		integral -= (hp/2)/sqrt(lambda + 2*intGamma);
      		next = intGamma - (Gamma((r+1)*hp)+Gamma((r+2)*hp))*hp/2;
      		integral -= (hp/2)/sqrt(lambda + 2*next);
      		for (int c=0; c<N; ++c)
        		(*initial)[r*N + c] = integral + t2;
    	}
    	int r = M-1;
    	intGamma -= (Gamma(r*hp)+ Gamma((r+1)*hp))*hp/2;
    	integral -= hp/sqrt(lambda + 2*intGamma);
    	for (int c=0; c<N; ++c)
      		(*initial)[r*N + c] = integral + t2;
   }
}
